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We present and demonstrate a method for optical homodyne tomography based on the inverse Radon trans- 
form. Different from the usual filtered back-projection algorithm, this method uses an appropriate polynomial 
series to expand the Wigner function and the marginal distribution and discretize Fourier space. We show that 
this technique solves most technical difficulties encountered with kernel deconvolution based methods and re- 
constructs overall better and smoother Wigner functions. We also give estimators of the reconstruction errors 
for both methods and show improvement in noise handling properties and resilience to statistical errors. 
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I. INTRODUCTION 

In quantum mechanics it is not possible to directly observe 
a quantum state \ip). In order to obtain full knowledge about 
|'0) it is necessary to accumulate measurement statistics of 
observables, such as position x or momentum p, on many dif- 
ferent bases. In quantum optics, this statistical measurement 
can be achieved by angle resolved homodyne measurement 
of the operator xg = x cos 9 + p sin 9 to acquire statistics of 
the squared modulus of the wave function |(a;e|-!/')p. Instead 
of the quantum state \tp), one is rather usually interested in 
reconstructing the more general density matrix p of the sys- 
tem. Fully equivalent to p, it is also possible to reconstruct 
the Wigner function W{q,p) from |(a;g|-0)p. However, the 
reconstruction of p or W{q,p) is not immediate and requires 
the reconstruction of the complex phase of the quantum sys- 
tem from the many angle resolved measurements. With the 
measurement of KxejV')!^, these two operations together are 
referred to as quantum homodyne tomography or optical ho- 
modyne tomography 1 1 ] . 

While some tomography algorithms reconstruct the for- 
mer density matrix, others rather reconstruct the latter Wigner 
function. Independently, tomography algorithms can be 
roughly classified into two species. Historically the first to be 
proposed and used for optical homodyne tomography, linear 
methods exploit and inverse the linear relationship between 
the experimentally measurable quantity Kxeli/;)! on one hand 
and p or W{q, p) on the other hand. Among them, the filtered 
back-projection algorithm IjJ, |2t] based on the inverse Radon 
transform |3] is the most commonly used. Similar in nature, 
there also exist methods based on quantum state sampling of 
individual components of the density matrix p with sample 
functions 14^,^. The linear methods, however, suffer in gen- 
eral from technical difficulties associated with the numerical 
deconvolution necessary to perform the linear inversion of the 
Radon transform (see Sec. II for details). In addition, they 
usually do not guarantee the physicality of the reconstructed 
state, the positivity of p. Finally they perform weakly against 
statistical noise and show numerical instabilities for higher 
frequency components and fine details of the reconstructed 
objects. Variational methods, such as the maximum entropy 
Og] and maximum likelihood [7J algorithms, were latter ap- 



plied to optical homodyne tomography to address these prob- 
lems. These methods can be designed to enforce the physi- 
cality of the reconstructed state and are usually more resilient 
to statistical errors. Since the reconstructed states are not de- 
fined constructively, an approximation procedure, typically it- 
erative, is used to achieve the reconstruction in practice [8]. 

Notice that in theory it is actually possible to bypass these 
numerical reconstructions and directly observe the Wigner 
function Wlq^p) with repeated measures of the parity oper- 
ator P — e*'^" where n is the number operator [9]. This mea- 
surement technique uses the link between the Wigner function 
value at point {q, p) and the expectation value of P for the dis- 
placed density matrix p 



W[q,p) = -tr D{-a)pD{a)e"^ 
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where Z) is the displacement operator and a = {q+ip)/^/2. A 
close tomography technique has been experimentally demon- 
strated in coupled systems of atoms and light |K|] . Unfortu- 
nately, a parity detector is a highly non-linear detector which 
can only be partially implemented for light beams with time- 
multiplexing and single photon detectors. Therefore with cur- 
rent state-of-the-art technologies in quantum optics, it is not 
possible to rely on count statistics alone for quantum state to- 
mography and one has to use optical homodyne tomography 
based on Gaussian measurements. 

While the linear methods look inferior to the variational 
methods, most of their associated problems are only technical 
in nature and can in principle be solved. In this paper we show 
that is it possible to use a linear reconstruction algorithm with 
better resilience to noise and better physical properties overall 
than the usual filtered back-projection method. The success 
of this approach lies in a systematic expansion of both the 
Wigner function W{q, p) and the marginal distribution p(x, 9) 
in polar coordinates. This circular harmonic expansion tech- 
nique has been applied in the past to other problems where 
the Radon transform plays a role in tomography 1 1 ll,fl2l], and 
here we adapt it to the quantum framework of optical homo- 
dyne tomography. In Sec. II we first review the basics of the 
inverse Radon transform and the usual filtered back-projection 
algorithms for optical homodyne tomography. In Sec. Ill we 
introduce the expansion method: we first conduct a spectral 
analysis of the angular components of p{x,9) and W{q,p); 
from this analysis we argue that a polynomial approximation 
is an efficient way to expand the radial components. In Sec. 



IV we give details about the implementation of the algorithm 
and also provide an estimator of the reconstruction errors. Us- 
ing our estimator we study the performances relatively to the 
filtered back-projection algorithm on simulated and experi- 
mental data sets. We complete this comparison with numer- 
ical studies of the distance between target and reconstructed 
quantum states. 



II. FILTERED BACK-PROJECTION 

In 1917, Radon introduces the integral transform 7?^ of two- 
dimensional functions integrated along straight lines and pro- 
vides the formula for the inverse transform 7?,"^ yfl- Today 
the Radon and inverse Radon transforms are ubiquitous in to- 
mography and find applications in many different area of sci- 
ence. The Radon transform is as well applicable to optical 
homodyne tomography. First we recall the definition of the 
observable operator xq of an homodyne measurement. 
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where Ug is the rotation operator in phase space, or phase- 
shifting operator The marginal distribution of the homodyne 
current p{x, 9) is then distributed according to the squared 
modulus of the wave function 



p{x,e) ^\{xe\ 



x\Ug\^){i,\U',\x), (2.2) 



where \xg) is the eigenvector of xg. The Radon transform TZ 
links the Wigner function W{q,p) of the quantum state \ip) 
and p{x, 9) the marginal distribution of the homodyne current 
with a projection of W{q, p) on a particular angle of observa- 
tion 9 In 



p{x,9) = n{w) 
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W{q,p)5{x — qcos9 — psin9)dqdp 



W{x cos 9 — p sin 9,xsm9 +p cos 9)dp. 

(2.3) 

In his original paper. Radon mathematically inverses his trans- 
form with the back-projection B of the derivative of the 
Hilbert transform H of p{x, 9) 

W{q,p) = ^B (-^n{p{x, 9)){y)}j , (2.4) 

where the back-projection operator S of a function f{x, 9) is 
the function F{q,p) defined by 

F{q,p)= f{qcoa9+psm9,9)d9. (2.5) 

Jo 

Expanding Eq. ( |2.4t we obtain the inversion formula 



{q cos 9 + p sin 9 ~ x)^ 



dxd9, 
(2.6) 



where V is the principal-value operator. Although exact, 
this expression is nevertheless unusable with experimen- 
tal data as the algebraic expression of p{x, 9) is unknown. 
However, the projection-slice theorem or Fourier slice theo- 
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Figure 1. Different transforms for different paths from p to W . 



rem lll4ll gives another reverse path from p{x, 9) to W{q,p) 
to work around the difficulties of the principal-value opera- 
tor (see FigdJ. \f p{k,9) and W{u,v) are, respectively, the 
one-dimensional and two-dimensional Fourier transforms of 
p(x, 9) and VF(q,p), then the projection-slice theorem states 
that 



p[k, 9) — W{k cos 6*, k sin t 



(2.7) 



Simply computing the Fourier transform p(fc, 9) from the 
measured data would seem like the most efficient way to ob- 
tain W{q,p) after a second inverse Fourier transform, but 
Eq. ( I2.7l i shows that it is necessary to interpolate W{u, v) 
in Fourier space, which leads to significant numerical difficul- 
ties 111511 ■ To avoid this interpolation Eq. ( 12.7b can be used to 
replace W{u, v) in the inverse Fourier transform of W{q,p) 
to obtain the inversion formula, 

W{q,p) = -- / p{x,9)K{qcos9+psm9-x)dxd9. 

(2.8) 
Here, the marginal distribution is convoluted with an integra- 
tion kernel K{x) and then back-projected into phase space, 
where K{x) is defined as the inverse Fourier transform of |fc| 
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(2.9) 



To use Eq. (I2.8l l in practice it is necessary to regularize K{x) 
and replace it with some numerical approximation. This is 
possible with the use of a window function g{k) such that the 
integral. 
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<7(fc)|fc|e*'=^dfc, 



(2.10) 



converges. The most common way to regularize Eq. ( 12.91 ) is to 
choose g{k) = l[_fc^ +fc^](fc) and introduce a hard frequency 
cutoff parameter kc so that 



K{x) 



1 



{cos(kcx) + kcX siiii{kcx) — 1) . (2.11) 
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Figure 2. Regularized integration kernel K{x) for different values of 



In practice, the choice of kc affects how much high fre- 
quency components of the Wigner function will get recon- 
structed. If kc is set too low the convolution in Eq .( 12.81 ) will 
filter out the fine physical details of the Wigner function. If fc^ 
is set too high, the convolution will introduce unphysical high 
frequency noise from the statistical errors in the measurement 
ofp{x, 9). Figure|2]shows the integration kernel K{x) for dif- 
ferent high frequency sensitivities. Choosing the right value of 
kc is a trade off between these two regimes. From Eqs. ( 12.7b 
and (12.81 ) it is also possible to insert other filter functions at 
different steps of the inversion to obtain modified algorithms 
with enhanced and more selective noise filtering properties. In 
any case the numerical implementation of Eq. ( I2.8l l will rely 
on deconvolution of the marginal distribution, an operation 
very sensitive to statistical noise. 



III. HARMONIC SERIES EXPANSION 

To numerically perform optical homodyne tomography, it is 
necessary at some point to apply an approximation procedure 
from the infinite dimensional space which features the un- 
known physical state to a finite dimensional space used to de- 
scribe the reconstructed state. In the filtered back-projection 
algorithm, the discretization is achieved by direct evaluation 
of W{xi,pi) on the set of points {{xi,pi)}i chosen to probe 
the phase space. Rather than this point-by-point reconstruc- 
tion, a discretization of another space should help to solve the 
numerical issues encountered in Sec. II. Since we are deal- 
ing with objects behaving like probability distributions, the 
statistical moments of p{x,d) and W{q,p) might be a solu- 
tion to the problem. In Ref.| 16], Ourjoumtsev et al. describes 
such a technique where they parametrize the Wigner function 
of a photon subtracted squeezed vacuum with the second and 
fourth moments of the marginal distribution p(x, 9). Gener- 
alizing this approach for any quantum state to higher order 
moments requires the use of the moment generating function 
<^e^^), where {x) is the expectation value of x with regards to 
p{x, 9). Superior to the moment generating function the char- 
acteristic function (e*^^) only needs the mean and variance 
to be defined to exist. This and the projection-slice theorem 
of Eq. ( 12.71 ) hint that Fourier space is a good candidate for an 
efficient discretization. 

We decompose our discretization procedure in two steps: 



(1) an angular harmonic decomposition with Fourier series; 

(2) a polynomial series expansion of the radial components. 
We express W{q, p) in radial coordinates (r, 0) and notice that 
M^(r, (/) + 27r) = W{r^ (j)). Therefore we write the radial part 
of M^(r, (/)) in terms of a Fourier series and we define the set of 
radial functions, or angular harmonic components {w„(r)}„ 
by 

Wn{r) = ^ / " W{r, <p)e-"'U<P, (3.1) 

which allows us to write W{r^ (j)), 

W{r,ci,)^ Y, ^nWe"^ (3.2) 

n— — oo 

with the symmetry relation w„ (r) — wl„ (r). The 2D Fourier 
transform W{u, v) of W{q, p) is written in radial coordinates, 

p + ca P+TT 

W{k,9)^ / T4^(r,0)e-*'''='=°''(^-*Vdr#, (3.3) 

Jo J-TT 

with the change of variables {u, v) -^ {k, 9). W{u, v) is re- 
lated to the Weyl function X(W)^') = tr(/3e-™'J+™P) by a 
simple 7r/2 rotation. 



W{u,v) = Xi-v,u), 
W{k,9) = Xik,9 + ^). 



(3.4) 
(3.5) 



We can easily write W in polar coordinates in terms of the 
angular harmonic components w„(r) of W{r, (p). 
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W{k,9)= Y^ / Wn{r)rdr 

n=-oo -^0 



With a Jacobi-Anger expansion of e*^ '^°^ "^ using Bessel func- 
tions J„, 



e'^''°^'f' = J2 «"^(^)' 



,m0 



(3.7) 



it is possible to conduct the angular integration in Eq. ( 13.6b to 
obtain the expression. 



oo />oo 

W{k,9)^2TT Y Hfe™'' / Wnir)Jn{kr 



)rdr. 



n— — oo 



(3.8) 
Notice that /„ Wn(r)Jn(kr)rdr is the n* order Hankel 
transform of Wn (r). 

In the same fashion, since p{x, 9 + 27r) = p{r, 9) we de- 
compose the marginal distribution as 



Pe 



{x)= Y c«(a;)e»^ 



(3.9) 



with the sets of radial functions c„ (x) defined by 



cn{x) = —l p{x,e)e-'^^'de. 



(3.10) 



Using the projection-slice theorem of Eq. (I2.7l i and the or- 
thogonality of e*"® on [— TT, +7r] we are able to write for every 
angular harmonic order n. 



■n /•+CXD ^oo 

-- / c„(x)e"'''^dx = / Wn{r)Jn{kr)rdr. 



(3.11) 



We have obtained a relation between, on one side the Fourier 
transform of the angular harmonics of p{x, 9), and on the 
other side, the Hankel transform of the angular harmonics of 
W{r, (p). If we inverse the Hankel transform with the orthog- 
onality relation, or closure relation of Bessel functions. 



1 



kdkJ„{kr)J„{kr') = -S{r - r') 



(3.12) 



we finally obtain 



w,.{r) = - 



Jn{kr)kdk 
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Cn{x)e-"'''dx. (3.13) 



At that point it would be natural to convey some radial de- 
composition of Wn{r) and Cn{x) . However, there is no simple 
way to achieve this. Looking at Eq. ( |3.13t , we notice that the 
Fourier transform of kJn{k), or at least J„(fc), should be in- 
volved in the process. The latter is written in terms of the 
Chebysheff 's polynomials of the first kind !"„ 

Uk)e-''^-dk = L_Z_ r„(x)l[„i,+i](x). 

-oo V ^ X 

(3.14) 
Equation (I3.14l i hints at the use of the polynomial series to 
achieve this radial decomposition. It is safe to assume for ap- 
plications that the Wigner function will only take nonzero val- 
ues from the origin up to a certain limit L > r. Since we are 
carrying the decomposition in polar coordinates what we are 
looking after is a polynomial family which is orthogonal on a 
disk of radius L. There are of course infinitely many such fam- 
ilies but one which proves to be particularly adequate to the 
task is the set of Zernike polynomials Z"(r. i^) — _R"(r)e""^ 
originally introduced for the study of o ptic al aberrations in 
lenses and other circular optical systems lUTIl . The polynomi- 
als are defined for s > \n\ > and s — |n| even. While the 
angular part gives straightforward orthogonality and fits with 
our previous approach using Fourier series, the radial compo- 
nents i?^" defined for t — \n\ > by 



{s-t)/2 

Rf-ir)^ J2 (-1)' 

fe=0 



(s - A:)! 



s-2/c 



fc!(^-fc)!(V-^)! ' 

(3.15) 
are orthogonal on [0, 1] with respect to the weight function r 
for all positive and negative orders n. 



R^{r)R^^,{r)rdr = —^< 
2(.s + 1) 



(3.16) 



Furthermore it turns out that the Radon transform of Zernike 
polynomials happens to have the simple expression, 

2 



n(R':{r)e"'^) = -v/l-x2C/,(x)e"^ 

s + 1 



(3.17) 



where Us(x) are the Chebysheff 's polynomials of the sec- 
ond kind lUSl \l% (see also the last paragraph of this section 
for a proof). The critical aspect for tomography lies in the 
fact that Us{x) is again an orthogonal polynomial family on 
[—1,1] with respect to the weight function y/1 — x^. In other 
words by finding a family of orthogonal polynomials whose 
Radon transform element by element is yet another family of 
orthogonal polynomials, we have in some sense diagonalized 
the Radon transform. The inverse Radon transform can also 
be exactly calculated and any technical difficulties associated 
with kernel functions or regularization immediately vanish. 

With the use of Eq. ( 13.161 ) we are eventually able to expand 
the angular harmonic functions w„ (r) on the n* order radial 
polynomials i?"(r). 



w;„(r)=^<i?,"(r). 



(3.18) 



s=0 



Given that R^{r) is non zero only when s > |n| > and s — 
|n| is even, we introduce the change of variable s -^ \n\ + 2m, 
re-index the sequence w^ and rewrite Eq. ( 13.18b 



^»(^) = Zl <'^|n|+2m('')- 



(3.19) 



m=0 



Putting Eqs. ( 13.191 ) and ( 13.2b together we obtain the complete 
expansion of W{r, </>) inside the unit disk D{Q, 1), 
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E <^|:|+2™('')e'"'- (3-20) 
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Notice fromEq. ( 13.151 ) that Rf^r) = -Rj"(r) which justi- 
fies the use of -R " , 2„j although w™ are in general complex 
constants. Applying the relation ( 13.171 ) on Eq. ( 13.20b . p(x, 9) 
is also written in terms of the coefficients w!" as 
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(3.21) 
To justify the use of Zernike polynomials and prove Eq. 
(13.17b . the relation. 
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RUr)Jn{rk)vdv = (_i)(™-«)/2:^ 
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(3.22) 



between Zernike polynomials and Bessel functions |i7| is es- 
sential. If we recall Eq. (13. lib , replace Wn (r) by its expansion 
on K'sir) in Eq. (13. 18b and cut the integration from +00 to 
unity, we obtain 
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-''"'^dx. 
(3.23) 



To finally obtain the complete inversion of TZ and the expan- 
sion of Cn{x) as in Eq. ( |3.21t . we only need to inverse the 
Fourier transform in Eq. ( 13.231 ) from the rhs to the Ihs and use 
the Fourier transform of Js{k)/k, 



Js+i{k) 



'dk = 



2r 



-C/,(a;)v/l-x2l[_i^+i](.T), 
(3.24) 



to obtain 
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(3.25) 
Notice that Eqs. ( |3.22t and ( |3.24| i close the link between 
Us{x) and _R",(r), the first two families of orthogonal func- 
tions used in the analysis, and the Bessel functions J„(/c) or- 
thogonal with respect to the weight function 1/fc, 



Js{k)Jt{k) 



dk 
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(3.26) 



In summary by identifying three families of orthogonal func- 
tions related together by the Radon transform TZ and the 
Fourier transform F, we have been able to find an expansion 
of the Wigner function W{q^ p) that allows to greatly simplify 
the technical difficulties of tomography with inverse Radon 
transform. 



IV. RECONSTRUCTION ALGORITHM 



A. The algorithm 




Figure 3. Comparison between polynomial series tomography (left 
panels: TV = 8, M — 30) and filtered back-projection tomography 
(right panels: fee = 8.0) for the state p = 0.8|1)(1| + 0.2|0)(0|. 
(a) J = 5 X 10^ (b) J = 20 X 10^ (c) J = 80 x 10^ (d) J = 
320 X 10^ (e) J = 5 X 10=*; (f) J = 20 x 10^ (g) J = 80 x 10^ 
(h) J — 320 X 10"^. All data sets have been synthetically generated 
with rejection sampling. 



The algorithm works in four steps: (1) choosing the size L 
of the reconstruction disk, (2) evaluating the coefficients w™, 
(3) choosing the cutoffs N and M of the angular and radial 
series, and (4) calculating W{r, (p). Step 1 is necessary for 
the orthogonal relations given in Sec. Hon [0, 1] and [—1,+!] 
to hold. In practice we have to normalize the marginal dis- 
tribution p{x, 9) — >■ p{x/L, 0)/L and the Wigner function 
W{r, (j)) -> W{r/L, (f>)/L. Step 2 is easily conducted by in- 
verting the relation (I3.21l i with the orthogonal Chebysheff's 
polynomials U\n\+2mix), 



w„ 



\n\ +2m + 1 
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dx U\n\+2m[X). 



The recurrence relation. 



Us+iix) = 2xUs{x) - Us^iix), 



(4.1) 



(4.2) 



allows one to efficiently calculate Us (x) for any s and any x 
given Uo{x) = 1 and Ui{x) = 2x. After obtaining the coef- 
ficients ui™ and choosing cutoff orders N and M, the Wigner 
function W{r, 0) is then approximated by the partial sums. 
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(4.3) 



Using the symmetry relation w™^ — (w™)*, we keep the real 
part of Eq. ( 14.31 ) and simplify the sum on n to 
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X (a™ cos(n0) + fe™ sin(n(/))) , 



(4.4) 



where we have defined w™ — (a™ + i6™)/2 for n > 1 and 
w™ — a™. Figures [3] and |4] show examples of reconstructed 
Wigner functions for a mixture of |0) and |1), and a thermal 
state respectively. In comparison to filtered back-projection 
tomography, polynomial series tomography converges faster 
with fewer numbers of experimental points J. The recon- 
structed Wigner functions also show less visible artifacts and 
are overall smoother. To evaluate efficiently -R™(r) we notice 
that W^(f) — rl"l and then use the recurrence relation ll20ll . 



-Rn+2(m+l)(^) 
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In contrast to setting the value of fc^, the values of N and 
M have a real physical meaning. This is a major advantage 




(a) 0.1 



Figure 4. Comparison between polynomial series tomography (left 
panels: N — 8, M — 30) and filtered back-projection tomography 
(right panels: kc — 8.0) for a thermal state of mean photon number 
(n) = 1. (a) J = 5 X 10^ (b) J = 20 x 10^ (c) J = 80 x 10^ (d) 
J = 320xl0^(e) J = 5xl0^(f) J = 20xl0^(g) J = 80x10^ 
(h) J = 320 X 10 . All data sets have been synthetically generated 
with rejection sampling. 




Figure 5. Effect of increased radial resolution on the stability 
of tomography of an experimentally measured photon subtracted 
squeezed vacuum (same data as in Ref. r27])- For all panels 
J = 1 X 10^. (a) Polynomial series tomography, A^ = 8, M = 20; 
(b) M = 30; (c) M = 40; (d) filtered back-projection tomography, 
kc = 7; (e) kc = 9; (f) fee = 11. 



of this method compared to the usual filtered back-projection 
algorithm. M will decide what will be the highest polynomial 
order of the radial features of W. Therefore it is equivalent 




Figure 6. Effect of A'^ and A/ on the convergence of polynomial 
series tomography. Same experimental data as in Fig. \5\ (a) Circular 
cut at constant r and effect of iV for A/ = 32, J = 2 x 10^. (b) 
Radial cut at constant (f> and effect of A/ for A^ = 10, J = 2 x 10^. 



to choosing the maximum photon number of the density ma- 
trix diagonal elements. N will set the resolution of the an- 
gular features of W, which decides how many off-diagonal 
components of the density matrix will be reconstructed. Fur- 
thermore it is easy to change N and AI after computing the 
coefficients w™. Figure |5] shows the effect of increasing M 
on the precision of polynomial series tomography. In compar- 
ison to filtered back-projection tomography when increasing 
the kernel sensitivity kc, increasing the radial resolution M 
does not produce artifacts in the Wigner function. Figure |6] 
further shows the effect of increasing N and M on the pre- 
cision of the tomography reconstruction of experimental data. 
While the angular components show quick convergence, the 
radial components require higher M values to be faithfully 
reconstructed. Figure Q illustrates the advantage of polyno- 
mial series tomography in radial resolution for quantum states 
with a higher number of photons. Both M and kc where set at 
values high enough to recover the original Schroedinger's cat 
state negativity at the origin of phase space. While the back- 
filtered projection shows numerical uinstability when kc is set 
high, the Wigner function reconstructed by polynomial series 
tomography is smoother at the equivalent resolution. 

Finally the value of Rn+2m ^nr = will be non-zero only 
for n — 0, therefore we have the useful formula to evaluate 
the Wigner function at the origin of phase space. 



M 



W^'(0,0)=5](-l) 



Vi, 



(4.6) 



which is similar to the formulation of W{0, 0) using the diag- 
onal elements of the density matrix. 



B. Unbiased error estimator 

To quantitatively compare our algorithm with the usual 
back-filtered tomography algorithm we give a consistent 
method to estimate the reconstruction error and obtain confi- 
dence intervals when calculating the value of W{q,p). If W 
and W" are the reconstructed value of W{q, p) with Eqs. ( |4.4| i 
and ( I2.8I 1 respectively, we call a^, and cr^„ the variance of 




Figure 7. Effect of increased radial resolution on the stability of to- 
mography of a Schroedinger's cat states with (n) — 3. For all panels 
J — 4: X 10*. (a) Original Wigner function; (b) polynomial series 
tomography, N — 8, M = 46; (c) filtered back-projection tomogra- 
phy, kc = 11. 



the reconstruction errors assuming they are distributed accord- 
ing to a Gaussian for both algorithms. We also assume that 
there are no systematic errors but only statistical eiTors. Let's 
assume an optical homodyne measurement set consists of J 
experimental points {{xj,9j)}j independently and identically 
distributed according to the underlying marginal distribution 
p{x,d). To begin with we give an estimator of <J^,,{q,p) 
for the usual filtered back-projection method using formula 
(I2.8I 1. To calculate the value of W at point {q,p), p{x, 9) will 
be replaced either by a binned histogram made from the data 
set {{xj , dj)}j, or by a sum of delta functions approximating 
pix,9) 

p{x, 6*) = ^ ^ (5(x - Xj) X d{9 - 9j). (4.7) 



In the latter case, the swap of p{x, 9) for expression (14.7b in 
Eq. (|2J]| leads to 



W"{q,p) 



1 ^ 

-—-y^K{qcos9,+psm9,-x,). (4.8) 



Since p{x,9) is a valid probability distribution W"{q,p) is 
nothing else than {K{q cos 9 + psin9 — x)) the expectation 
value of the kernel function. Therefore Eq. J4.8I ) can be re- 
garded as a Monte Carlo integral where the expectation value 
of the kernel function is calculated by randomly sampling K 
according to the distribution p{x, 9). In other words, the op- 
tical homodyne tomography with filtered back-projection is 
in effect an analogical Monte Carlo integration where the ho- 
modyne measurement plays the part of the random number 
generator In that familiar case the statistical properties of the 
reconstruction error are well known. First of all we are as- 
sured of the unbiased convergence of the sum in Eq. ( 14.81 ). 
The central limit theorem also states that the eiTor will indeed 
converge to a Gaussian distribution of zero mean and whose 
standard deviation aw" {q, p) for J experimental points is 



aw"{q,p) = ctk/VJ-I, (4.9) 

which exhibits a 1/a/J rate of convergence, and where gk = 
y/{K^) — {K)'^/2n. By using the approximations. 



1 ^ 

(K) « — 2^ K{q cos 9j + p sin 9j — Xj), 



(4.10) 



3=0 
J 



{K"^) w -^K'^{qcos9j + piiin9j - Xj), (4.11) 

we can actually estimate ax in a straightforward way easy to 
include in the implementation of Eq. ( 14.8b . 

The same analysis for the coefficients {w™} yields the re- 
construction sum. 



2m 



27r2 



1 ^ 



• 1+2' 



,(x,/L)e-^"^Vi- (4.12) 



As previously errors are Gaussian distributed for every coeffi- 
cient ui™ with a l/\/j rate of convergence. If a quantity Y is 
calculated through the measure of the variables {yi]i<i with 
the formula. 



Y = fiyi,---,yi), 

then the variance a-y of Y can be approximated by 



(4.13) 



(4.14) 
where a"^ ~ (xy) — {x){y). Using Eq. (14.41 ) we can ap- 
ply this formula to estimate the variance aw anywhere in 
phase space, but because of its simple formulation thanks to 
Eq. ( 14.61 ). we will only study it at the origin (0, 0); 
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(4.15) 
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Figure 8. Estimation of (tw{Q, 0) with filtered back-projection to- 
mography (plain line) and polynomial series tomography (dotted 
lines). (a)p = 0.8|1)(1| +0.2|0){0|; (b) thermal state with (n) = 1; 
(c) photon subtracted squeezed vacuum (same data as in Fig[5j. 



Notice that in this case the variance estimator formula of Eq. 
(14.14b is not an approximation anymore due to the linear com- 
bination nature of Eqs. (14.4b or (14.6b . We can compute an esti- 
mate of CTa'" when computing the coefficients w™ in the same 
way we did with Eqs. (14.10b and ( 14.11b . Figure |8] shows es- 
timation of the reconstruction errors for different states using 
Eq. ( 14.91 ) and ( 14. 15b . We have found that the value of k^ has 
very little influence on aw" at the center of phase space. On 



the contrary M has a strong influence on aw' (0, 0). However, 
as was shown in Figsl3]|4]and|5] far from the origin the poly- 
nomial series tomography algorithm shows less uncertainties. 
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Figure 9. Effect of M on the convergence of W {Q, 0) and the mag- 
nitude of a"iy/(0, 0). (a) Thermal state with (n) = 1, rejection sam- 
pling, (b) Experimental photon subtracted squeezed vacuum state 
(same data as in Fig|5j. 

We also assumed the convergence error due to finite trun- 
cation N and M of the expansion to be smaller than the sta- 
tistical error itself. This can be checked in the algorithm by 
iteratively calculating a^, (0, 0) for increasing values of M 
and stop when the magnitude of the A/* and last coefficient 
Wq^ is less than u^, (0, 0) (see Fig. |9]l. This technique can be 
repeated independently for every point of phase space {q,p), 
and different values of N and M can even be used for different 
points of phase space. 



C. Monte Carlo error estimation 

Independently from the estimators of the previous para- 
graph, we also use Monte Carlo simulations to generate many 
synthetic data sets and evaluate the reconstruction errors. This 
method is easily applied if we know precisely which state \ip) 
is under investigation. For example, we can choose a known 
density matrix or Wigner function and calculate the associated 
marginal distribution p{x, 6). From this marginal distribution 

we generate K synthetic data sets of J points {{xj^Oj)}:^ 
using, for example, rejection sampling. With the algorithm 
of our choice we repeat the tomography reconstruction and 
calculate a set of K Wigner function {VF'^'''-'}fc. Finally for a 



(a) 



simulated 
estimated 




100 
K 



(b) 



simulated 
estimated 




(a) 1 r 




(b) 1 



1000 10000 

J 




10000 100000 

J 



Figure 10. Comparison between Monte Carlo simulation and direct 
estimation of crvv(0, 0). Black curves are the estimation of aw (0,0) 
with Monte-Carlo simulation using K data sets. Dashed curves are 
the direct estimation of aw{0,0) using Eqs. l |4.9l l and l l4.15t for 
the K'*' data set. (a) Data sets of J = 10^ points generated using 
rejection sampling for the state 0.8il)(l| + 0.2|0)(0i. (b) Data sets 
of J = 10^ points generated with bootstrapping resampling from 
the same experimental data as Fig. |5] (i) Filtered back-projection 
tomography with kc = 7; (ii) polynomial series tomography with 
M = 10; (iii) M = 20; (iv) M = 30; (v) M = 40. 



given point of phase space {xo,po), we calculate the average 
value Wo of the set {W'^'^^k- 

1 ^ 

^o = 1fE^^''^(^"'Po)' ^4.16) 



k=l 



and obtain an estimate of the error a-gy at point (xo,po) by 

, K 2 

^VF = 17E(^''^(^0'Po)-^o) • (4.17) 



k=l 



Since it is a Monte Carlo based simulation, every quantity 
shows again a 1/VK convergence rate. 

With experimental data, we can sample p{x, 9) only once 
and therefore we need a technique to generate the synthetic 
data sets after the experimental measurement. Resampling is 
the easiest approach and here we estimate the reconstruction 
error of experimental data sets with the bootstrapping resam- 
pling method ||22I1 . The results of both techniques are illus- 
trated in Fig. [To] and overall there is a good agreement be- 
tween the estimated values of Monte Carlo simulations and 
the predicted value of aw{0, 0) using Eq. (14. 9t or (14.15b . 



D. Distance to a target state 

To conclude this comparative study of polynomial series 
expansion and filtered back-projection-based tomography, we 
numerically estimate in this final paragraph the distance be- 
tween some original target quantum state and reconstructed 
states using both algorithms. For this purpose we will consider 
one distance for the Wigner function and one distance for the 
density matrix. We use the L2 Euclidian distance dL2{-, •) for 



Figure 11. Estimation of the distance between the target thermal state 
of mean photon number (n) = 1 and reconstructed quantum states 
averaged over 1000 samples of J data points for different tomogra- 
phy settings, (a) L2 distance (dL2(VFtarge(, Womo))- (b) Frobenius 
distance (dF(/5tiugct,/5tomo)). 
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Figure 12. Estimation of the distance between the target state 
0.8|l)(lj + 0.2i0)(0| and reconstructed quantum states averaged 
over 1000 samples of J data points for different tomography set- 
tings, (a) L2 distance {dL2{Wtmget, VKtomo)). (b) Frobenius distance 

[dp (Ptarget , Ptomo ) ) . 
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Figure 13. Estimation of the distance between the target odd 
Schroedinger's cat state oc |a) — | — q) with (h) = 3 and 
reconstructed quantum states averaged over 1000 samples of J 
data points for different tomography settings, (a) L2 distance 

(dL2(Wtargei, VKtomo)). (b) Frobenius distance (dF(Aaigei,Aomo)). 



the Wigner function defined by 



dL2 {WA,WB)=(f J dxdp \Wa{x,p)~ Wb 



{x,p)\' 



1/2 



(4.18) 
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and with the Frobenius norm \\.\\f defined by 



1/2 



\Ay = V^A^= El^ 



■y I ' 



(4.19) 



we define a distance di?(., .) for density matrix as 

dFipA,PB) = \\pA-pB\\. (4.20) 

First we choose a target state and derive its exact Wigner func- 
tion Wtaiget and density matrix pt^ixst- We then evaluate the dis- 
tances from the target state according to Eqs. ( 14.181 1 and ( 14.20b 
using as before Monte Carlo sampling techniques. Rather than 
averaging a reconstructed state over many simulated data sets, 
we average the distance computed over many reconstructed 
states and estimate the numbers: 

(dL2(Wtarget,Wtomo)) and (c^F (Ptaiget, Aomo)) • (4.21) 



Numerical simulation results are shown in Figs. [TTIfT3] for. 
respectively, a thermal state with (n) = 1, a mixture of vac- 
uum and one-photon state 0.8|1)(1| + 0.2|0)(0|, and an odd 
Schroedinger's cat state with (n) — 3. In agreement with 
the previous results on tomography uncertainties, we observe 
that polynomial series expansion tomography performs bet- 
ter than filtered back-projection for these two first cases. In 
the case of the Schroedinger's cat state ex \a) — \ — a), both 
distances behave differently for higher J and tend to reach a 
precision limit which depends on the tomography algorithm 
and settings. Although the exact cause of this saturation is 
unknown, we believe it is due to the significantly more com- 
plex structure of the Schroedinger's cat state. According to 
our simulations, it seems to depend only on the radial and an- 
gular precision settings, more precisely on parameters AI, N, 
and kc- In this case again, polynomial series expansion proves 
to reach a higher precision level than filtered back-projection 
for a relevant range of tomography settings. To conclude this 
paragraph, it is interesting to notice that in the case of the 
dL2{-, ■) distance there is an intrinsic limitation on the preci- 
sion of polynomial series expansion tomography due to the 
circular geometry of the reconstruction space |23]. This could 
be the reason for the saturation phenomenon visible in Fig. 

m 



V. CONCLUSION 

We have shown and demonstrated a technique for optical 
homodyne tomography based on polynomial series expansion 



of the Wigner function. In Sec. II we have given the basis 
of the usual filtered back-projection algorithm and explained 
the main reason for its weak performances against statisti- 
cal noise. We have also introduced the projection-slice the- 
orem and the relation between phase space, Fourier space and 
the marginal distribution. In Sec.III we have shown that it 
is possible to link three families of orthogonal functions be- 
tween these three spaces to decompose p{x, 9) the marginal 
distribution, W{q,p) the Wigner function, and their Fourier 
transforms. We have shown that the Radon transform pre- 
serves the orthogonality of these families and therefore takes 
an especially simple form in this case. In Sec. IV we have 
explained and applied to experimental and simulated data the 
most straightforward implementation of that technique with a 
direct linear estimation of the coefficients of the polynomial 
series expansion. We have also provided estimators of the re- 
construction errors and shown that it performs better than fil- 
tered back-projection tomography with respect to reconstruc- 
tion artifacts and statistical errors. More precisely, polynomial 
series tomography is superior with fewer experimental data 
points and when higher radial resolution is needed for higher 
photon number states. These results are confirmed when look- 
ing at the distance between a chosen target state and states re- 
constructed with both tomography techniques. Furthermore 
this technique exploits the projection slice theorem directly 
and therefore is faster than convolution based filtered back- 
projection. Finally we remark that it is in principle possible to 
use the maximum likelihood technique to find the set of coef- 
ficients ?«J^ that maximizes the probability of measuring the 
experimentally measured data set. 
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